#drafted by Jeremy VanDerWal ( jjvanderwal@gmail.com ... www.jjvanderwal.com )
#GNU General Public License .. feel free to use / distribute ... no warranties

################################################################################
###get the command line arguements
args=(commandArgs(TRUE)); for(i in 1:length(args)) { eval(parse(text=args[[i]])) } #evaluate the arguments
#sample data
# es="RCP85"
# gcm="ukmo-hadgem1"
################################################################################
###define a couple functions
#define a function to convert data & lat/lon to ascii grid
matrix2ascii = function(lat,lon,tdata) {
	lats = unique(lat); lats = sort(lats)
	longs = unique(lon); longs = sort(longs)
	cellsize = min(c(diff(lats), diff(longs)))
	nc = ceiling((max(lats) - min(lats))/cellsize) + 1
	nr = ceiling((max(longs) - min(longs))/cellsize) + 1
	out.asc = as.asc(matrix(NA, nr = nr, nc = nc), xll = min(longs), yll = min(lats), cellsize = cellsize)
	return(put.data(cbind(lon,lat,tdata), out.asc))
}
#define a funtion to remove NA data by buffering with averages
buff.na = function(tt.asc) {
	dims = dim(tt.asc) #get the dimensions
	tmat = matrix(NA,nrow=dims[1]+2,ncol=dims[2]+2) #create a temporary matrix
	tmat[2:(dims[1]+1),2:(dims[2]+1)] = tt.asc[,] #copy the data
	out = tt.asc #make a copy of the data
	while(!all(is.finite(out))) { #cycle through until all NA values are filled
		for (yy in 2:(dims[1]+1)) { #cycle through rows
			for (xx in 2:(dims[2]+1)) { #cycle through columns
				if (!is.finite(tmat[yy,xx])) { #if it is an NA data
					tmat[yy,xx] = mean(as.vector(tmat[(yy-1):(yy+1),(xx-1):(xx+1)]),na.rm=TRUE) #get the mean of surrounding cells
					out[yy-1,xx-1] = tmat[yy,xx] #copy the mean to the output ascii
				}
			}
		}
	}
	return(out)
}
################################################################################
###start doing the work
library(climates); library(SDMTools) #define the libraries needed

wd = '~/Climate/CIAS/Australia/5km/';setwd(wd) #define and set the working directory
fd = '~/Climate/CIAS/global.20120322/' #define the locations of the future data

#read in the current observed data at 5km
pos.obs = read.csv('baseline.76to05/base.positions.csv',as.is=TRUE) #read in the positions for the observed data
baseasc.obs = read.asc.gz('baseline.76to05/base.asc.gz') #read in the base asc for the observed data
prec.obs = as.matrix(read.csv('baseline.76to05/monthly.pr.csv',as.is=TRUE))
tmin.obs = as.matrix(read.csv('baseline.76to05/monthly.tmin.csv',as.is=TRUE))
tmean.obs = as.matrix(read.csv('baseline.76to05/monthly.tmean.csv',as.is=TRUE))
tmax.obs = as.matrix(read.csv('baseline.76to05/monthly.tmax.csv',as.is=TRUE))

#read in current climate at 50km
pos.cur = read.csv(paste(fd,'base.positions.csv',sep=''),as.is=TRUE) #read in the positions for the modelled data
baseasc.cur = read.asc.gz(paste(fd,'base.asc.gz',sep='')) #read in the base asc for the current modelled data
prec.cur = as.matrix(read.csv(paste(fd,'monthly_csv/current.1976.2005/pre.matrix.csv',sep=''),as.is=TRUE))
tmin.cur = as.matrix(read.csv(paste(fd,'monthly_csv/current.1976.2005/tmn.matrix.csv',sep=''),as.is=TRUE))
tmean.cur = as.matrix(read.csv(paste(fd,'monthly_csv/current.1976.2005/tmp.matrix.csv',sep=''),as.is=TRUE))
tmax.cur = as.matrix(read.csv(paste(fd,'monthly_csv/current.1976.2005/tmx.matrix.csv',sep=''),as.is=TRUE))

#read in teh future data given the ES & GCM of interest
prec.fut = as.matrix(read.csv(paste(fd,'monthly_csv/',es,'/',gcm,'/pre.matrix.csv',sep=''),as.is=TRUE))
tmin.fut = as.matrix(read.csv(paste(fd,'monthly_csv/',es,'/',gcm,'/tmn.matrix.csv',sep=''),as.is=TRUE))
tmean.fut = as.matrix(read.csv(paste(fd,'monthly_csv/',es,'/',gcm,'/tmp.matrix.csv',sep=''),as.is=TRUE))
tmax.fut = as.matrix(read.csv(paste(fd,'monthly_csv/',es,'/',gcm,'/tmx.matrix.csv',sep=''),as.is=TRUE))

years = unique(as.numeric(substr(gsub('pre','',colnames(prec.fut)[-c(1:2)]),1,4))) #get the years in the datasets

#calculate the anomolies
prec.anom = prec.fut; tmin.anom = tmin.fut; tmean.anom = tmean.fut; tmax.anom = tmax.fut #make a copy of the data
for (year in years) { 
	cois = grep(year,colnames(prec.anom)) #define the columns to be working with 
	prec.anom[,cois] = prec.fut[,cois] / (prec.cur[,-c(1:2)]+1) #calculate anomoly as a proportion
	tmin.anom[,cois] = tmin.fut[,cois] - tmin.cur[,-c(1:2)] #calculate anomoly as absolute difference
	tmean.anom[,cois] = tmean.fut[,cois] - tmean.cur[,-c(1:2)] #calculate anomoly as absolute difference
	tmax.anom[,cois] = tmax.fut[,cois] - tmax.cur[,-c(1:2)] #calculate anomoly as absolute difference
}

#downscale and apply the anomolies
out.prec = out.tmin = out.tmean = out.tmax = cbind(pos.obs[,c('lat','lon')],matrix(NA,nrow=nrow(pos.obs),ncol=ncol(prec.anom)-2)) #setup the output data matrices
colnames(out.prec) = colnames(prec.anom); colnames(out.tmin) = colnames(tmin.anom); colnames(out.tmean) = colnames(tmean.anom); colnames(out.tmax) = colnames(tmax.anom); #define the column names
rois = which(pos.cur[,'lat'] <= -6 & pos.cur[,'lat'] >= -52 & pos.cur[,'lon'] >= 108 & pos.cur[,'lon'] <= 160) #define the rows of the anomolies that are needed for interpolation
xout = getXYcoords(baseasc.obs)$x; yout = getXYcoords(baseasc.obs)$y #define the output x & y values
for (ii in 3:ncol(prec.anom)) { cat(colnames(prec.anom)[ii],'..',ii,'\n') #cycle through each of the columns
	##precip
	tmat = buff.na(matrix2ascii(c(prec.anom[rois,'lat'],-51.75),c(prec.anom[rois,'lon'],159.75),c(prec.anom[rois,ii],NA))) #get the matrix around australia and fill NA values
	tt = interp2grid(tmat,xout,yout) #do the interpretation
	out.prec[,ii] = tt[cbind(pos.obs[,'row'],pos.obs[,'col'])] * (prec.obs[,(((ii-3)%%12)+1)+2]+1) #apply the anomoly
	##tmin
	tmat = buff.na(matrix2ascii(c(tmin.anom[rois,'lat'],-51.75),c(tmin.anom[rois,'lon'],159.75),c(tmin.anom[rois,ii],NA))) #get the matrix around australia and fill NA values
	tt = interp2grid(tmat,xout,yout) #do the interpretation
	out.tmin[,ii] = tt[cbind(pos.obs[,'row'],pos.obs[,'col'])] + tmin.obs[,(((ii-3)%%12)+1)+2] #apply the anomoly
	##tmean
	tmat = buff.na(matrix2ascii(c(tmean.anom[rois,'lat'],-51.75),c(tmean.anom[rois,'lon'],159.75),c(tmean.anom[rois,ii],NA))) #get the matrix around australia and fill NA values
	tt = interp2grid(tmat,xout,yout) #do the interpretation
	out.tmean[,ii] = tt[cbind(pos.obs[,'row'],pos.obs[,'col'])] + tmean.obs[,(((ii-3)%%12)+1)+2] #apply the anomoly
	##tmax
	tmat = buff.na(matrix2ascii(c(tmax.anom[rois,'lat'],-51.75),c(tmax.anom[rois,'lon'],159.75),c(tmax.anom[rois,ii],NA))) #get the matrix around australia and fill NA values
	tt = interp2grid(tmat,xout,yout) #do the interpretation
	out.tmax[,ii] = tt[cbind(pos.obs[,'row'],pos.obs[,'col'])] + tmax.obs[,(((ii-3)%%12)+1)+2] #apply the anomoly	
}
csvdir = paste(wd,'monthly_csv/',es,'/',gcm,'/',sep=''); dir.create(csvdir,recursive=TRUE) #create the output csv directory
write.csv(out.prec,paste(csvdir,'pre.matrix.csv',sep=''),row.names=FALSE) #write out the data
write.csv(out.tmin,paste(csvdir,'tmn.matrix.csv',sep=''),row.names=FALSE) #write out the data
write.csv(out.tmean,paste(csvdir,'tmp.matrix.csv',sep=''),row.names=FALSE) #write out the data
write.csv(out.tmax,paste(csvdir,'tmx.matrix.csv',sep=''),row.names=FALSE) #write out the data

#create the bioclim surfaces / mxe files
for (year in years) { cat(year,'\n')
	##get the subset climate data for the time of interest
	tpre = as.matrix(out.prec[,grep(year,colnames(out.prec))])
	ttmin = as.matrix(out.tmin[,grep(year,colnames(out.tmin))])
	ttmean = as.matrix(out.tmean[,grep(year,colnames(out.tmean))])
	ttmax = as.matrix(out.tmax[,grep(year,colnames(out.tmax))])
	
	out = bioclim(ttmin,ttmax,tpre,ttmean) #create the bioclim variables
	colnames(out) = paste('bioclim_',sprintf('%02i',as.numeric(gsub('bioclim_','',colnames(out)))),sep='') #change the coloumn names
	#out = cbind(out.prec[,1:2],out) #append location info
	
	ascdir = paste(wd,'bioclim_asc/',es,'_',gcm,'_',year,'/',sep=''); dir.create(ascdir,recursive=TRUE) #define and create output directory
	for (ii in colnames(out)) {
		tasc = baseasc.obs; tasc[cbind(pos.obs[,'row'],pos.obs[,'col'])] = out[,ii] #put the data in ascii format
		write.asc(tasc,paste(ascdir,ii,sep=''))
	}
	mxedir = gsub('bioclim_asc','bioclim_mxe',ascdir); dir.create(mxedir,recursive=TRUE) #define and create directories
	system(paste('java -mx1024m -cp ~/maxent.jar density.Convert ',ascdir,' asc ',mxedir,' mxe',sep='')) #convert to mxe
	system(paste('gzip ',ascdir,'*.asc',sep='')) #gzip the data
};
